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We study the organization of topological defects in a system of nematogens confined to the two- 
dimensional sphere {S'^). We first perform Monte Carlo simulations of a fluid system of hard rods 
(spherocylinders) living in the tangent plane of S'^ . The sphere is adiabatically compressed until 
we reach a jammed nematic state with maximum packing density. The nematic state exhibits four 
+ 1/2 disclinations arrayed on a great circle rather than at the vertices of a regular tetrahedron. This 
arises from the high elastic anisotropy of the system in which splay (i^i) is far softer than bending 
(K3). We also introduce and study a lattice nematic model on S'^ with tunable elastic constants 
and map out the preferred defect locations as a function of elastic anisotropy. We establish the 
existence of a one-parameter family of degenerate ground states in the extreme splay-dominated 
limit KzjKx 00. Thus the global defect geometry is controllable by tuning the relative splay to 
bend modulus. 
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Ordered arrays of microscopic structures on curved 
interfaces provide a promising route to fabricating 
nanoscale or mesoscale building blocks (mesoatoms) that 
may in turn form molecules and bulk materials. One class 
of mesoatoms is provided by particles self-assembling on 
spherical droplets in liquid-liquid emulsions. The parti- 
cles may be isotropic or shaped. Ordered structures on 
spherical interfaces always possess topological defects [1] . 
The location and detailed arrangement of such defects 
are both important since defects are distinctive regions 
which may be functionalized to create directional bonds 
akin to atomic bonds. This has been nicely illustrated 
recently by the synthesis of divalent gold nanoparticles 
coated with self-assembled stripes of phase-separated lig- 
ands whose two polar defects can be functionalized. The 
divalent nanoparticles subsequently link spontaneously 
to form chains [2]. 

Spherical nematics have a local 2-fold inversion sym- 
metry and elementary disclinations of both half-integer 
and integer strength. Since the total disclination strength 
on the 2-sphere is two O |4] , nematic ground states may 
possess four +1/2 defects, two +1 defects or two +l/2s 
and one +1. Functionalization of the defects in the first 
case could lead to tetravalent mesoatoms with sp^-like di- 
rectional bonding [3 . The structure and arrangement of 
defects in such a thin nematic shell with variable thick- 
ness has recently been studied both experimentally [5 
and theoretically [6 . The analyses so far have, however, 
been limited to the one Frank constant approximation, in 
which the bending stiffness K-^ and the splay constant K\ 
are equal [7 . Although fluctuations always drive these 
elastic moduli to the same value in the long wavelength 
limit, we are necessarily working at finite volume on the 
compact 2-sphere. It is essential, therefore, to explore the 
effect of differing bend and splay moduli on the structure 
of defects in the ground state. 

To this end we have performed Monte Carlo simula- 
tions of hard rod fluids confined to the tangent plane 



of the 2-sphere. Adiabatically shrinking the sphere in- 
creases the packing density and leads to a jammed state 
with nematic order and four +1/2 disclinations. Since 
like-sign defects repel, one might expect them to be max- 
imally separated at the vertices of a regular tetrahe- 
dron [3l H]. We find instead that the four defects lie 
approximately on a great circle. This can be understood 
as arising from the high bending stiffness compared 
to the splay stiffness constant K\. We also analyze a 
coarse-grained model of a spherical nematic with tun- 
able Frank constants and map out the global pattern of 
defects as a function of the anisotropy K^/Ki. In the 
limit K2,/Ki ^ 00, we show that the system exhibits a 
one-dimensional continuum of degenerate ground states 
in which the four defects form a rectangle of arbitrary 
aspect ratio. 

Our Monte Carlo simulation is performed in the 
isobaric-isothermal (constant-NPT) ensemble with rod- 
like particles interacting via hard core repulsion We 
use spherocylinders with length L and diameter D (as- 
pect ratio L/D). Isotropic initial states of N rods are 
prepared by centering each randomly oriented rod on a 
node of a spherical mesh of a 2-sphere of initial radius Ri. 
We perform ~ 10^ Monte Carlo cycles, each of which con- 
sists of a translational and orientational trial move for all 
N rods together with a volume- compression move. The 
compression rate must be low enough to avoid prema- 
ture jamming at low density. The simulation stops when 
the system is jammed and thus cannot be further com- 
pressed. 

We simulated a system of iV = 1082 rods with as- 
pect ratio L/D = 15. The final radius reached was 
Rf/D = 37.15 starting from Ri/D = 100. The final 
states are found to be nematic with four +1/2 discli- 
nations. This contrasts with the corresponding infinite 
system in 2D flat space where the densest packing state 
is crystalline [9j. At these system sizes, therefore, the 
spatial curvature of the surface frustrates crystalline or- 
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FIG. 1: (Color online) The ground-state configuration for 
spherical nematic ordering of 1082 rods with L/D = 15 is 
shown as the central lines of rods for clarity, together with 
maps of the local nematic order parameter. The two views 
are: (a) near defects and (b) nematic bulk. The radius of the 
compressed sphere is Rf / D = 37.15. The measurement of 
the order parameter along the great circle connecting defects 
is plotted as a function of angular distance in radians in (c). 
The exact locations of the four defects lying near one great 
circle are indicated in (d). 



dering. Crystalline order is nevertheless expected to set 
in for larger system sizes, where the radius of curvature 
becomes much larger than the rod length. In simulations 
with short rods {L/D < 7) we find, in contrast, smectic- 
like domain structures at maximum density. 

To characterize the nematic order we measure the lo- 
cal order parameter tensor Qap^ a symmetric, traceless 
second- rank tensor, by locally averaging the orientation 
of m{V) rods n{a) within a small domain V (typically 
with an area of order L^): 
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The nematic order parameter S is defined to be twice 
the positive eigenvalue of Qa^- Maps of S are shown in 
Figs. [TJa) and (b). The local order parameter along the 
great circle connecting the defects is also plotted versus 
angular distance in Fig. [l}^c). The maximum value of S 
is about 0.93, in contrast with flat space where S reaches 
1 [9]. The reduced maximum on the sphere results from 
the frustration due to the nonzero curvature enclosed in 
the region V. The defect cores are identified as the loci of 



minima of the order parameter. Four +1/2 disclinations 
are clearly observed to lie on a single great circle, to 
within ±0.15rad, as illustrated in Fig. [TJd). 

The planarity of the defect array can be explained by 
the strong elastic anisotropy. It is known that, for hard 
rod systems, the bend constant K^, diverges as the rod 
density increases while the splay constant Ki is almost 
independent of density [TOl [TTj . This leads to ^ Ki 
at high packing density. To check this we determine the 
Frank constants by fitting the local director field sur- 
rounding a defect to a numerical minimization of the 
Frank free energy [l2]. We find that Ks/Ki > 20. 

Consider first planar nematics in the strong anisotropy 
limit. It is known that +1 disclinations admit both pure 
bending texture or pure splay texture J2| ^3 • When the 
bending and splay constants are different one of these two 
textures will correspond to a minimum of the energy and 
the other a maximum. The energy of a +1 disclination 
is therefore given by F^+i) - 7rmin(i^i,i^3)ln(i^/ao), 
where R is the system size and ao is the defect core 
size. The energy of a +1/2 defect, by contrast, depends 
on both Ki and Ks. In the one Frank constant limit 
Ki = Ks = i^, the energy of a +1/2 disclination is one- 
quarter that of a +1 disclination: ^ \n{R/ao). 
A +1 disclination will therefore unbind to a pair of 
+1/2 disclinations which may then separate. In the 
splay-dominated case Ki <C i^s, in contrast, the energy 
of +1/2 disclination is half that of a +1 disclination: 
- ^\n{R/ao). Splitting a +1 defect therefore 
gains no energy and the total energy of two +1/2 discli- 
nations is independent of their separation. All bending 
modes are frozen for a 2D nematic in flat space with 
Ks = 00. It can be shown in this limit that the in- 
tegral curve of the director field is a geodesic (straight 
line). This implies that the 2D nematic director texture 
is completely determined by its configuration on a one- 
dimensional curve: there is no bulk degree of freedom 
that one can vary. In Fig. [2|a) we illustrate the forma- 
tion and separation of a pair of +1/2 pure-splay defects 
by global surgery on a +1 pure- splay defect. The +1 






(a) 



FIG. 2: (Color online) (a) In the limit of infinite K3/K1, the 
separation of a +1/2 disclination pair in the plane costs no 
energy, (b) On the 2-sphere, an infinite number of states with 
four +1/2 disclinations on a great circle can be generated from 
a state with two +1 disclinations by cut-and-rotate surgery. 
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defect is pulled apart horizontally and the intermediate 
region filled with uniform nematic texture. Clearly it is 
also possible to pull apart the defect pair vertically with- 
out free energy cost. By contrast, any infinitesimal local 
deformation of the state in Fig. [2|a) generates bending 
deformations and is therefore forbidden energetically. 

We now extend this argument to nematic order on the 
2-sphere with infinite bending stiffness. The constraint of 
no bending again translates into the requirement that the 
integral curves of the director follow geodesic lines (great 
circles). Consider an initial state for which the director 
field follows lines of longitude everywhere. This state, 
with one +1 disclination at each pole, is clearly bending- 
free and is therefore a ground state in the limit Ks oo. 
Any local deformation of this state introduces bend and 
is therefore forbidden. There are, however, global manip- 
ulations that cost zero energy. As illustrated in Fig.[2|b), 
we can cut the sphere into two hemispheres along a great 
circle that contains both +1 disclinations, and rotate one 
hemisphere by an arbitrary angle a. This surgery divides 
each +1 disclination into two +1/2 disclinations and the 
resultant four +1/2 disclinations all lie on the great circle 
along which we cut. Note that the director field is ev- 
erywhere smooth up to first order derivatives on the cut 
after the surgery except at the disclination cores. 
The elastic free energy of the post-surgery state is there- 
fore independent of the rotation angle a. We have thus 
identified a continuous manifold of degenerate low en- 
ergy states, parameterized by the angle a, for spherical 
nematics in the pure- splay limit Ks oo. The energy 
difference between the four +1/2 disclination state and 
the original state with two +1 disclinations can only come 
from the defect core energies and these do not grow with 
radius R. It is thus negligible for large system sizes. 

Starting from a pure splay state and rotating the direc- 
tor at every point by 7r/2 in the tangent plane, we obtain 
a pure bend state. It can also be shown that the splay 
and bending constants Ki and Ks interchange under this 
rotation. A system with K^/Ki <C 1 should, therefore, 
also exhibit a similar one-parameter family of degenerate 
ground states with a co-planar configuration of defects. 
This is verified below. We note that the same cut-and- 
rotate surgery has been applied to spherical smectics (for 
which Ks/Ki is also ^ 1) by Chantawansri et al; [16^ and 
Blanc and Kleman [17]. 

To understand better how the global configuration of 
defects depends on the elastic anisotropy Ks/Ki, we next 
analyze a lattice nematic model on S'^ in which the two 
Frank constants can be continuously varied. We first 
construct a triangular lattice on the 2-sphere with N = 
762 vertices. As a result of the Gauss-Bonnet theorem, 
twelve of these vertices have coordination number five 
and form a regular icosahedron. It turns out, however, 
that these crystalline disclinations do not significantly 
affect the location of the nematic defects ^8J. On every 
vertex a of this lattice we define a nematic director n^, 
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FIG. 3: (Color online) The effect of the anisotropy J3/J2 on 
the configuration of defects. As the anisotropy is turned on 
the defects shift from a tetrahedral geometry to a coplanar 
great circle alignment. The angular deviation of the defect 
positions from the great circle are measured as a function of 
the anisotropy. 



which is constrained to the local tangent plane. The 
Hamiltonian is given by 

H = J2Y. Tr(A^a - N^){Na - N^) (2) 

(a,6) 

+ J3 ^ {eab'Na-eacmNa-N^)-{Na-N,), 

(a,6/c) 

where Na is the dyadic tensor fiafia^ while e^b is the unit 
vector pointing from vertex a to vertex h. In Eq. ([2| (a, 6) 
means summation over all pairs of nearest-neighbor sites 
h and c that are both nearest neighbors of site a. 

The lattice model Eq. ([2| is manifestly invariant un- 
der spatial inversion of each director n^- The two body 
term with coefficient J2 is precisely the 2D Maier-Saupe 
lattice model of nematics, which is known to reduce, at 
large scales, to the Frank free energy with equal elastic 
constants Ki = Ks. The three body term with coeffi- 
cient J3 couples the directors explicitly to the underly- 
ing lattice and renders Ki and Ks different. It can be 
shown [19 that at large scales the lattice model Eq. ([2| 
reduces to the following Frank free energy for spherical 
nematics: 
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where D is the covariant derivative and t is the unit tan- 
gent vector perpendicular to the director h. We note that 
in going from Eq. ([2| to Eq. (|3|, however, terms of higher 
order (in Dn) have been dropped. These terms, though 
irrelevant in the renormalization group sense, may be 
quantitatively important for small system sizes. The pa- 
rameter Js in Eq. ([2| turns out to be proportional to 
{Ks — Ki)/{Ks + 2Ki). The coefficient of proportional- 
ity, however, depends on microscopic details [19 . 
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FIG. 4: (Color online) At J3/J2 = 1.2, the energy, as a func- 
tion of the rotation angle a, is flat within ±2%. 



again at zero temperature. As shown in Fig. |4j the en- 
ergy of all low energy states thus obtained, for different 
angular separation a, agree to within ±2%. This small 
difference is likely due to finite size effects. 

In this letter we have studied the global configuration 
of defects in the ground state of a spherical nematic using 
two different models. We show that the tetrahedral con- 
figuration crosses over to a great- circle configuration as 
we increase the anisotropy of the elastic constants. Our 
work therefore explicitly demonstrates that defect posi- 
tions can be controlled by varying the elastic anisotropy. 
This result should be relevant to designing and fabricat- 
ing mesoscopic molecules and bulk materials by attach- 
ing ligands to functionalized defect sites. The work of 
MJB and HS was supported by the NSF through Grants 
DMR-0219292 and DMR-0305407. The work of XX was 
supported by the American Chemical Society through 
Grant PRF 44689-G7. MJB thanks Robin Selinger for 
discussions. 



Starting from a random state, we determine the ground 
state of the model Eq. ([2| using simulated annealing. We 
use the same method to identify the nematic defects as 
we used for the hard rods simulation. We first simulate 
the one Frank constant case for which J3 = and find 
that the four defects form a regular tetrahedron as ex- 
pected [31 [4] . A non-zero three-body coupling J3 leads to 
many metastable states. To avoid them, we start from 
the minimum energy configuration for the one Frank con- 
stant case and turn on the parameter J3 at sufficiently 
low temperature. We then cool down again to reach the 
ground state. We find that the nematic texture becomes 
increasingly splay (bend) dominated as J3 increases (de- 
creases) from zero. We also measure the angular devia- 
tions of the defect locations from the circumscribed great 
circle, determined by a least squares fit. As shown in 
Fig. [3j the defects gradually shift toward a great circle as 
IJ3I increases. When J3/J2 > 1.2 or < — 0.6, the direc- 
tor field develops a local instability, from which we infer 
that the splay or bend constant (at short scales) becomes 
negative. Nevertheless, the defect array is not strictly 
planar at J3/J2 ^ 1-2 or — 0.6. We suspect that this is 
due to the higher order terms that have been dropped in 
going from the discrete model Eq. ([2| to the continuous 
model Eq. (|3| , which in principle can renormalize the two 
Frank constants and make them scale dependent. Such 
effects, however, are expected to diminish at large system 
sizes. 

The lattice nematic model can be used to display ex- 
plicitly the ID manifold of degenerate ground states. For 
J3/J2 = 1-2, we start with a low energy state obtained by 
simulation of Eq. Q in which the four defects lie approx- 
imately on the equator with an angular separation of 90°. 
We then force two defects to move within the great cir- 
cle by rotating all the directors in the lower hemisphere. 
Starting from this new state we then minimize Eq. (H 
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